subroutine gridearnings
use params
use mini
implicit none

!variables to get stationary distribution on epsilon
integer, parameter :: itmax_eps=10000
integer :: iter_eps
double precision, parameter :: tol_eps=1e-7
double precision :: d_eps, put_eps0(neps)

double precision, dimension(nz,nz) :: pi_z
double precision, dimension(nz) :: grid_z

double precision, dimension(neta,neta) :: pi_eta
double precision, dimension(neta) :: grid_eta, put_eta0

double precision, dimension(nnu,nnu) :: pi_nu
double precision, dimension(nnu) :: grid_nu

double precision, dimension(nx,nx) :: pi_x

integer :: id
integer :: jeta, jz, jnu, idp
double precision :: temp_u1(nu1,nu1)
    
integer, parameter :: n_earnings=neps*nu1
integer :: i_earnings, id_eps(n_earnings), id_u1(n_earnings), iter_simulate
double precision, dimension(n_earnings,neps,nu1) :: omega_earnings, omega_earnings_p

double precision :: earnings, earnings_p, mass_p, earnings_initial
double precision :: avg_log_earnings_change_1, var_log_earnings_change_1
double precision :: avg_log_earnings_change_5, var_log_earnings_change_5

double precision :: avg_log_earnings_perm, var_log_earnings_perm
double precision :: dummy

double precision :: avg_log_earnings, var_log_earnings, std_log_earnings
double precision :: avg_log_earnings_p, var_log_earnings_p, std_log_earnings_p
double precision :: avg_log_earnings_both

integer :: cnt_log_earnings_p
double precision, allocatable, dimension(:,:) :: dist_log_earnings_p
double precision :: log_earnings_p_90, log_earnings_p_10
integer :: find

double precision, dimension(6) :: avg_log_earnings_t, var_log_earnings_t, std_log_earnings_t, avg_log_earnings_both_t

! if ( rank .eq. 0 ) then

    if ( nu1 .eq. 1 ) then
    
    grid_u1=1
    pi_u1=1
    
    else
    
        call tauchen (nu1,rho_z,mu_z,var_z,grid_u1,temp_u1,1*one,3*one)
        
    grid_u1=exp(grid_u1)
    pi_u1=temp_u1(1,:)
!     grid_u1=grid_u1/sum(grid_u1*pi_u1)

    end if
    
    if ( rank .eq. 0 ) then
    
        open(unit=1,file='output_grid_u1.txt')            
        write(1,*)'u1 pi_u1'
    
            do iu1=1,nu1
            write(1,*)grid_u1(iu1),pi_u1(iu1)
            end do
        
        close(1)
        
    end if

    if ( nz .eq. 1 ) then

    pi_z=1
    grid_z=1

    else
    
        call tauchen (nz,rho_z,mu_z,var_z,grid_z,pi_z,one,1*one)
    grid_z=exp(grid_z)

    end if

    if ( neta .eq. 1 ) then
    
    pi_eta=1
    grid_eta=1
    
    else
    
        call tauchen (neta,rho_eta,mu_eta,var_eta,grid_eta,pi_eta,1.*one,4.*one)
    grid_eta=exp(grid_eta)

!         call tauchen (neta,rho_eta,mu_eta,var_eta,grid_eta,pi_eta,one,13.*one)
!         call TransitionOneStep(mu1,sigma1,mu2,sigma2,zero,zero,p1,p2,rho_eta,grid_eta,pi_eta)    
!         
!     grid_eta=exp(grid_eta)
    
    end if
    
    if ( nnu .eq. 1 ) then
    
    pi_nu=one
    grid_nu=1
    
    else
    
        call tauchen (nnu,rho_nu,mu_nu,var_nu,grid_nu,pi_nu,one,2.*one)
    grid_nu=exp(grid_nu)        
    
    end if
    
pi_x=one

put_eps0=zero

    do inu=1,nnu    
    put_eps0((inu-1)*nz*neta*nx+1)=pi_nu(inu,inu)
    end do

pi_nu=zero

    do inu=1,nnu
    pi_nu(inu,inu)=one
    end do
    
    id=0
    
    do inu=1,nnu
    do iz=1,nz
    do ieta=1,neta
    do ix=1,nx   
    
    id=id+1
    id_nu(id)=inu
    ix_ieps(id)=ix    
    grid_eps(id)=grid_nu(inu)*grid_eta(ieta)*grid_z(iz)
    
    idp=0
    
        do jnu=1,nnu
        do jz=1,nz    
        do jeta=1,neta
        do jx=1,nx                        
        idp=idp+1
        pi_eps(id,idp)=pi_nu(inu,jnu)*pi_eta(ieta,jeta)*pi_z(iz,jz)*pi_x(ix,jx)
        end do
        end do
        end do
        end do
        
    end do
    end do    
    end do
    end do
        
d_eps=1
iter_eps=0
put_eta0=one/neta

    do while ( iter_eps .le. itmax_eps .and. d_eps .ge. tol_eps )

    put_eps=matmul(put_eps0,pi_eps)
    put_eta=matmul(put_eta0,pi_eta)
    d_eps=max(maxval(abs(put_eps0-put_eps)),maxval(abs(put_eta0-put_eta)))
    put_eps0=put_eps
    put_eta0=put_eta

    iter_eps=iter_eps+1

    end do

    if ( rank .eq. 0 .and. iter_eps .ge. itmax_eps ) then
    print*,'put_eps did not converge ',d_eps
!         call exit
    end if

! avg_eps=sum(put_eps(:)*grid_eps(:))    
avg_eps=zero

    do iu1=1,nu1
    do ieps=1,neps
    avg_eps=avg_eps+pi_u1(iu1)*put_eps(ieps)*grid_u1(iu1)*grid_eps(ieps)
    end do
    end do

! grid_u1=grid_u1/avg_eps
! grid_eps=grid_eps/avg_eps   

grid_beta(1)=one
grid_beta(2)=.5
put_beta(1)=one-quasi_hyperbolic_share
put_beta(2)=quasi_hyperbolic_share
    
   if ( rank .eq. 0 ) then
       
        open (unit=1,file='output_grid_eps.txt')
        write(1,*)'eps put_eps'
            do ieps=1,neps
            write(1,'(300F14.10)')grid_eps(ieps),put_eps(ieps)
            end do
        close (1)
        
        open (unit=1,file='output_grid_u1_eps.txt')
        write(1,*)'w*u1*eps u1*eps mass'
            do iu1=1,nu1
            do ieps=1,neps
            write(1,'(300F14.10)')w*grid_u1(iu1)*grid_eps(ieps)/avg_eps,grid_u1(iu1)*grid_eps(ieps)/avg_eps,pi_u1(iu1)*put_eps(ieps)
            end do
            end do
        close (1)
        
        open (unit=1,file='output_put_eta.txt')
        write(1,*)'put_eta'
            do ieta=1,neta
            write(1,'(300F14.10)')put_eta(ieta)
            end do
        close (1)
                
   end if    
   
!count for log permanent earnings distribution
! omega_earnings=zero
! i_earnings=0
! 
!     do ieps=1,neps
!     do iu1=1,nu1
!     i_earnings=i_earnings+1
!     omega_earnings(i_earnings,ieps,iu1)=put_eps(ieps)*pi_u1(iu1)    
!     end do
!     end do   
! 
! cnt_log_earnings_p=0
!     
!     do iter_simulate=1,2
!     
!     omega_earnings_p=zero
!         
!         do i_earnings=1,n_earnings
!         do ieps=1,neps
!         do iu1=1,nu1
!         
!             if ( omega_earnings(i_earnings,ieps,iu1) .gt. zero ) then
! 
!                 do jeps=1,neps
!                 do ju1=1,nu1
! 
!                 mass_p=omega_earnings(i_earnings,ieps,iu1)*pi_eps(ieps,jeps)*pi_u1(ju1)                        
!                 
!                     if ( mass_p .gt. zero ) then
!                 
!                         if ( iter_simulate .eq. 2 ) then
!                         cnt_log_earnings_p=cnt_log_earnings_p+1        
!                         end if
!                                 
!                     omega_earnings_p(i_earnings,jeps,ju1)=omega_earnings_p(i_earnings,jeps,ju1)+mass_p
!                 
!                     end if
!                 
!                 end do
!                 end do
! 
!             end if
!             
!         end do
!         end do 
!         end do
! 
!     omega_earnings=omega_earnings_p
!     
!     end do

!     allocate (dist_log_earnings_p(cnt_log_earnings_p,2))

!initial distribution
omega_earnings=zero
i_earnings=0

    do ieps=1,neps
    do iu1=1,nu1
    i_earnings=i_earnings+1
    id_eps(i_earnings)=ieps
    id_u1(i_earnings)=iu1
    omega_earnings(i_earnings,ieps,iu1)=put_eps(ieps)*pi_u1(iu1)    
    end do
    end do
    
avg_log_earnings_change_1=zero
var_log_earnings_change_1=zero

avg_log_earnings_change_5=zero
var_log_earnings_change_5=zero

avg_log_earnings_perm=zero
var_log_earnings_perm=zero

avg_log_earnings=zero
var_log_earnings=zero

avg_log_earnings_p=zero
var_log_earnings_p=zero

avg_log_earnings_both=zero

! dist_log_earnings_p=zero
! cnt_log_earnings_p=0

avg_log_earnings_t=zero
var_log_earnings_t=zero
avg_log_earnings_both_t=zero

    do iter_simulate=1,6
    
    omega_earnings_p=zero
        
        do i_earnings=1,n_earnings
        do ieps=1,neps
        do iu1=1,nu1
        
            if ( omega_earnings(i_earnings,ieps,iu1) .gt. zero ) then

            earnings_initial=w*grid_eps(id_eps(i_earnings))*grid_u1(id_u1(i_earnings))/avg_eps
            earnings=w*grid_eps(ieps)*grid_u1(iu1)/avg_eps                
            
                do jeps=1,neps
                do ju1=1,nu1

                mass_p=omega_earnings(i_earnings,ieps,iu1)*pi_eps(ieps,jeps)*pi_u1(ju1)                        
                
                    if ( mass_p .gt. zero ) then
                
                    earnings_p=w*grid_eps(jeps)*grid_u1(ju1)/avg_eps
            
                        if ( iter_simulate .eq. 1 ) then                
                    
                        avg_log_earnings_change_1=avg_log_earnings_change_1+log(earnings_p/earnings)*mass_p
                        var_log_earnings_change_1=var_log_earnings_change_1+((log(earnings_p/earnings))**2.)*mass_p   
                                            
                        avg_log_earnings=avg_log_earnings+log(earnings)*mass_p
                        var_log_earnings=var_log_earnings+(log(earnings)**2)*mass_p

                        avg_log_earnings_p=avg_log_earnings_p+log(earnings_p)*mass_p
                        var_log_earnings_p=var_log_earnings_p+(log(earnings_p)**2)*mass_p

                        avg_log_earnings_both=avg_log_earnings_both+(log(earnings)*log(earnings_p))*mass_p

                        end if
                
                        if ( iter_simulate .eq. 5 ) then
                        avg_log_earnings_change_5=avg_log_earnings_change_5+log(earnings_p/earnings_initial)*mass_p
                        var_log_earnings_change_5=var_log_earnings_change_5+((log(earnings_p/earnings_initial))**2.)*mass_p                
                        end if
                
                        if ( iter_simulate .eq. 2 ) then
                        dummy=(log(earnings_initial)+log(earnings)+log(earnings_p))/3.                
                        avg_log_earnings_perm=avg_log_earnings_perm+dummy*mass_p
                        var_log_earnings_perm=var_log_earnings_perm+(dummy**2)*mass_p                                
!                         cnt_log_earnings_p=cnt_log_earnings_p+1        
!                         dist_log_earnings_p(cnt_log_earnings_p,1)=dummy
!                         dist_log_earnings_p(cnt_log_earnings_p,2)=mass_p                        
                        end if
                                
                    avg_log_earnings_t(iter_simulate)=avg_log_earnings_t(iter_simulate)+log(earnings)*mass_p
                    var_log_earnings_t(iter_simulate)=var_log_earnings_t(iter_simulate)+(log(earnings)**2)*mass_p
                    avg_log_earnings_both_t(iter_simulate)=avg_log_earnings_both_t(iter_simulate)+(log(earnings_initial)*log(earnings_p))*mass_p

                    omega_earnings_p(i_earnings,jeps,ju1)=omega_earnings_p(i_earnings,jeps,ju1)+mass_p
                
                    end if
                
                end do
                end do

            end if
            
        end do
        end do 
        end do

        if ( iter_simulate .eq. 1 ) then
        std_log_earnings_change_1=sqrt(var_log_earnings_change_1-avg_log_earnings_change_1**2)
        std_log_earnings=sqrt(var_log_earnings-avg_log_earnings**2)
        std_log_earnings_p=sqrt(var_log_earnings_p-avg_log_earnings_p**2)
        correlation_lag_1=(avg_log_earnings_both-avg_log_earnings*avg_log_earnings_p)/(std_log_earnings*std_log_earnings_p)
        end if

        if ( iter_simulate .eq. 5 ) then
        std_log_earnings_change_5=sqrt(var_log_earnings_change_5-avg_log_earnings_change_5**2)
        end if

        if ( iter_simulate .eq. 2 ) then
        std_log_earnings_perm=sqrt(var_log_earnings_perm-avg_log_earnings_perm**2)
        end if
        
    std_log_earnings_t(iter_simulate)=sqrt(var_log_earnings_t(iter_simulate)-avg_log_earnings_t(iter_simulate)**2)

    omega_earnings=omega_earnings_p
    
    end do

    do iter_simulate=1,5
    correlation_lag_t(iter_simulate)=(avg_log_earnings_both_t(iter_simulate)-avg_log_earnings_t(iter_simulate)*avg_log_earnings_t(iter_simulate+1))&
    /(std_log_earnings_t(iter_simulate)*std_log_earnings_t(iter_simulate+1))
    end do
    
! dist_log_earnings_p(:,2)=dist_log_earnings_p(:,2)/sum(dist_log_earnings_p(:,2))
! 
!     call sort(dist_log_earnings_p,cnt_log_earnings_p)
! 
! find=0
! 
!     do iter_simulate=1,cnt_log_earnings_p
!     
!         if ( find .eq. 0 .and. sum(dist_log_earnings_p(1:iter_simulate,2)) .ge. .1*one ) then
!         find=1
!         log_earnings_p_10=dist_log_earnings_p(iter_simulate,1)
!         end if
!         
!     end do
!     
! find=0
! 
!     do iter_simulate=1,cnt_log_earnings_p
!     
!         if ( find .eq. 0 .and. sum(dist_log_earnings_p(1:iter_simulate,2)) .ge. .9*one ) then
!         find=1
!         log_earnings_p_90=dist_log_earnings_p(iter_simulate,1)
!         end if
!         
!     end do
!     
!     deallocate (dist_log_earnings_p)
! 
! log_earnings_p_90_10=log_earnings_p_90-log_earnings_p_10

! end if

end subroutine gridearnings
